Critical Temperature Curve in the BEC-BCS Crossover 
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The strongly-correlated regime of the BCS-BEC crossover can be realized by diluting a system of 
two-component fermions with a short-range attractive interaction. We investigate this system via 
a novel continuous-space-time diagrammatic determinant Monte Carlo method and determine the 
universal curve T c /sf for the transition temperature between the normal and the superfiuid states 
as a function of the scattering length with the maximum on the BEC side. At unitarity, we confirm 
that T c /e F = 0.152(7). 

PACS numbers: 



In the area of ultracold gases, the problem of the 
crossover between the Bardeen-Cooper-Schrieffer pair- 
ing and the Bose-Einstein condensation of composite 
molecules (the so-called BCS-BEC crossover) has re- 
cently received a lot of theoretical and experimental at- 
tention A dilute two-component Fermi gas, where the 
inter-particle distance is much larger than the interaction 
range, features a remarkable universality at low temper- 
atures. Since the interaction is completely described by 
the s-wave scattering length a, the only physically rel- 
evant coupling parameter is k — 1/kpa, where kp is 
the Fermi momentum. One thus obtains a unified and 
universal description of systems as diverse as ultracold 
fermionic gases in magnetic or optical traps [l[ , fermions 
in optical lattices, inner crusts of neutron stars 0, 
and, plausibly, excitonic condensates 

In the limit k — > — oo, the Fermi gas is described by 
the BCS theory, while for k — > +oo the fermions pair into 
compact bosonic molecules which then form a BEC state 
below the critical temperature. Separating these extreme 
states is a strongly correlated regime which features the 
so-called unitary point k = 0. At unitarity, the scattering 
length is infinite and the interaction thus drops out of 
the relations between different thermodynamic potentials 
making these relations formally identical to those of a 
non-interacting Fermi gas [f|. On the experimental side, 
using the technique of a (wide) Feshbach resonance in a 
system of cold atoms, one can traverse the whole range 
of parameter k from the BEC to the BCS limit 

Despite considerable recent investigation, the quanti- 
tative description of the BEC-BCS crossover is far from 
being complete, even for the simplest case of the equal 
mixture of two components. Due to the strongly corre- 
lated nature of the problem, analytical mean-field-type 
calculations (e.g. [1,0, Hi) unavoidably involve approxi- 
mations, the accuracy of which is difficult to access unless 
the exact result is known. Renormalization group treat- 
ments can be carried out as expansions in either e = 4 — d 
Q , or 1 /Nf (where Nf is the number of fermion species) 



(ToL [lH , but the applicability of these calculations to the 
physically relevant case of d — 3 and Nf — 2 is not 
known a priori. 

Numerical studies of fermionic systems are computa- 
tionally demanding and further complicated by the need 
to study the limit of small densities to access the uni- 
versal regime. Some numerical techniques avoid the 
fermionic sign problem with a help of uncontrollable ap- 
proximations. The restricted path-integral Monte Carlo 
(R-PIMC) [lH relies on a variational ansatz for the nodes 
of the density-matrix. In the dynamical mean-field the- 
ory (DMFT) approach of Ref. [l3| the physics of ex- 
tended paired states is reduced to that of a single site cou- 
pled to the self-consistently defined environment. Fortu- 
nately, the unpolarized Fermi gas with contact attrac- 
tion is an exceptional case which can be addressed by 
sign-problem-free determinant methods without uncon- 
trollable systematic errors 14, 151, Moreover, the de- 
terminant diagrammatic MC approach (DDMC) for lat- 
tice fermions [HI is completely free of any systematic 
error. In simulations of the negative-?/ Hubbard model 



14j with an appropriate extrapolation to zero filling at 



the unitary point we previoulsy obtained accurate results 
for the critical temperature of the superfiuid transition, 
T c /sf = 0.152(7), in the units of the Fermi energy ef- 
This result, however, did not agree with the estimate ob- 
tained by Ref. [16[ from a visual inspection of the caloric 
curve, using the standard auxiliary field approach |17|. 

To efficiently study the critical temperature curve away 
from the unitarity point and to verify that the final re- 
sults are model independent thereby also resolving the 
controversy on T c /ef at unitarity - we develop a DDMC 
technique for continuous space and time. We can now ef- 
ficiently simulate models with a simple parabolic disper- 
sion relation and have completely eliminated lattice cor- 
rections. In this Letter, we first discuss the new scheme 
and how to obtain an independent systematic-error-free 
value for T c /ef at unitarity. We are able to reach densi- 
ties almost 20 times smaller than those typically accessi- 



2 



T/s. 

C r 








0.25 




i 


A.----A- A a 


0.20 






BEC 


0.15 








0.10 








0.05 








0.00 


^BCS 

i 


i 


i 


-1 





1 a: 



scale Zq such that 



FIG. 1: (Color online.) The universal results for the critical 
temperature in the units of the Fermi energy plotted versus 
k — \ jkpa (circles). The solid lines for negative and positive 
K represent the limiting behavior of the BCS theory (with 
the Gorkov-Melik-Barkhudarov correction) and the ideal BEC 
respectively. For reference, we also plot non-universal results 
for hard-sphere (triangles) and soft-sphere (squares) bosons 
Ref. El. 



ble with the auxiliary field determinant method [16], This 
allows us to perform a reliable extrapolation to the uni- 
versal limit yielding T c /ef = 0.152(7), in perfect agree- 
ment with our previous value Next, we explore the 
critical temperature at finite values of l/kpa. Our re- 
sults, shown in Fig. [TJ fix the general shape of the uni- 
versal curve T c /ep versus l/kpa. The main feature is a 
substantial maximum of T c /ep on the BEC side of the 
crossover. 

Our specific model is described by the Hamiltonian 

- U J dx^x^Jto^x^Cx) , (1) 



<r=U' 



where 4 r (T (x) is the fermion field operator (er =f, |), x is a 
continuous three-dimensional coordinate, fi is the chem- 
ical potential, U < is the contact interaction strength, 
and K is the kinetic energy operator, Ke 1 ^ = £ke lkx , 
with £k being the single-particle dispersion. 

The scattering length a is given by the sum of the 
vacuum ladder diagrams [l8| leading to (h = 1) 



m 
Ana 



= U' 



dk 1 

(2^)3 2^ ' 



(2) 



where m is the fermion mass. For the continuous space 
model with — k 2 /2m an ultraviolet regularization of 
Eq. J2J is required. Keeping in mind comparison with 
Ref. [16j | , where the parabolic dispersion with an ultravi- 
olet cutoff was used, we introduce a microscopic length 



k 2 /2m , k < 2ir/l 
oo , k > 2tt/Iq 



yielding 



m/ina = U~ 1 - U; 



= —nlo/m 



(3) 



(4) 



It is straightforward to generalize the DDMC method 
for resonant fermions P3] to the continuous model JTJ). 
One starts by expanding the partition function Z = 
Tre~" H , where j3 = 1/kgT, in powers of U. The re- 
sulting Fcynman diagrams consist of four-point inter- 
action vertices connected by free single-particle prop- 
agators Go? . A diagram of a given order p is de- 
scribed by the space-time configuration of the vertices 
S p = {(x^Tj), j = l,...,p)} (r S [Q,0\ is the imagi- 
nary time) and the topology of propagator lines connect- 
ing them without integration over the vertex positions — 
the latter is done by the Monte Carlo sampling process. 
Next, one observes [l9| that the sum over all topologies 
is given by det A T (S p ) det (S p ), where A CT is the px p 

matrix, Afj(S p ) = G^fc — Xj,Tj — Tj). In the case of 
equal densities of the spin components, the weight of a 
configuration S p is positive definite: 



olP(p,S p ) = (-[/y |detA(Sp) 



; II dT 3 dx J 



(5) 



The partition function Z — X^Lo Is ^ IS calculated 
stochastically according to the standard Metropolis- 
Rosenbluth 2 -Teller 2 algorithm ensuring that configura- 
tions S p are generated with the probability density given 
by Eg. ([5]). The Monte Carlo updates are based on a 
worm-algorithm for the four-point correlation function 
[H G 2 (x,r;x',r') = <T T P(x, r)Pt(x', t')> , where T r in- 
dicates time-ordering, P(x, t) = ^(x, r)^ , ;(x } r) is the 
pair annihilation operator, and (■ • • ) is the thermal aver- 
age. The asymptotic value of ffdrdr' Ga(x, r;x' ,t') as 
x — x'| — » co is proportional to the condensate density. 

Up to statistical errors, the DDMC calculations yield 
exact results for a finite system - in our case a cubic box 
of a linear size L with periodic boundary conditions. An 
efficient way of finding T c in the thermodynamic limit 
L — > oo is to employ the technique of Binder crossings 
Q for R = L 1+ " / dxdx'drdr'G 2 (x,r;x',r')/(/3L 3 ) 2 
(where r] w 0.038 for the 3D U(l) universality class), 
as discussed in detail in Ref. [la ]. It is expected that at 
the critical point R becomes scale invariant. By analyz- 
ing the crossings of the family of R(L, 0) curves one can 
obtain T c with an accuracy of a fraction of percent with a 
relatively small number of particles. The thermodynamic 
limit of the number density is obtained from a linear ex- 
trapolation of n(L) as a function of l/L. An example of 
the finite-size analysis for a typical set of parameters is 
given in Fig. [21 
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FIG. 2: (Color online.) Finite-size analysis for c = 0.83, 
fj, — 0.36, corresponding to U ~ —7.519 and a « 1.52 (in the 
units of m = 1/2, l = 1) yielding /3 C = 1.290(8). The error 
bars are one standard deviation and were calculated using the 
blocking method. Inset: the thermodynamic limit value of the 
number density is obtained via a linear fit of n(L) vs. 1/L. 
In this case, nil = 0.119(2), which results in £ = 0.492(3), 
TP /e ( p = 0.335(6), and k = 0.432(3). 

In order to obtain the universal answer for T c /ef one 
finally has to take the limit of ( = rt 1 / 3 ^ — > by extrap- 
olating numerical data for tP/sP to the dilute limit. 
We keep lo constant and take the limit by lowering the 
chemical potential /j, and diluting the system. One can 
show [l5| that the leading-order corrections should be lin- 
ear in C < 1: TP /e ( p = T c /e F + const x(+ o(C). The 
calculation strategy is as follows: at unitarity (k = 0), 
we fix U = U* according to Eg. (0| and perform a se- 
ries of simulations for different values of fi, yielding a 
set of tP/sP. Then, the universal value of the criti- 
cal temperature follows from the linear extrapolation of 
T { P/eP to C -> 0. 

To obtain T c /e F away from the resonance, the pro- 
cedure has to be modified. Taking the dilute limit for 
each value of k ^ requires that a — » oo in such a 
way that 1/kpa tends to a fixed finite value k. We 
note that the universal value of the chemical potential 
obeys (j,(T c )/ef = 2mg(n), with some function g(n), or, 
equivalently, lim^o (t^) (T c )a 2 = g{n)/K 2 implying that 
for each k one has to keep fia 2 — const. Substituting 
a 2 = fi/c into (HJ) gives 

where the upper (lower) sign corresponds to the BEC side 
a > (BCS side a < 0). We thus pick a value of c and 
perform a series of simulations for smaller and smaller 
values of // with U from Eq. Each simulation yields 
a finite-^ estimate for the critical temperature tP(c), 



density n^(c) and k^(c). After linear extrapolations 
to C — > we determine the physical value of T c /sf and 
the corresponding value of k. 
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FIG. 3: (Color online.) The extrapolation of the simulation 
results to the universal limit £ = n l ' 3 lo 0. The proce- 
dure yields T c /e F = 0.152(9), 0.202(9), and 0.252(15) for 
k = l/k F a = (squares), 0.217(2) (circles), and 0.474(8) 
(triangles) correspondingly. For comparison, we also plot our 
results for the Hubbard model (open squares) adapted from 
Ref. [13 • The estimate of Ref. [lo ] at k = (obtained for fi- 
nite ( ps 0.93) is shown by the diamond. Solid lines are linear 
fits. 

In Fig. [3] we show results for the critical temperature as 
a function of £. For comparison and consistency analysis 
of T c at unitarity, we also plot the data for the Hubbard 
model [3] a s a function of the filling factor v, which 
plays the same role as C m the present model. Note that 
the non-universal corrections to T c /e F in ( turn out to 
be positive and much smaller (at unitarity) than for the 
Hubbard model. The former fact is important for the 
simulation efficiency, since the computational complexity 
of the DDMC technique scales as (/3UN) 3 , where N is the 
number of fermions, and it is advantageous to simulate 
at higher temperatures. 

It is important to note that at high densities the 
TP/sP curves are almost constant and the true asymp- 
totic low-C behavior develops only below ( w 0.75. For 
a reliable extrapolation it is crucial to vary the density 
by at least an order of magnitude, and we did so by di- 
luting the system down to n an 0.04/Zq (£ rj 0.35), where 
we were limited by the low values of the absolute criti- 
cal temperature itself. Unfortunately, no dilute-limit ex- 
trapolation was performed in Ref. [1 61 ] (their value for 
C w 0.93 is shown by the diamond in Fig. [3]). The total 
simulation time required to obtain this set of data was 
approximately 10 6 CPU hours on Opteron-class worksta- 
tions. 
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At unitarity, the extrapolated result for tJ^/e^P yields 
an answer which is in perfect agreement with T c /ef — 
0.152(7) obtained independently from the v — » extrap- 
olation of the Hubbard model data [3]. In the latter 
case, the universal value is approached from below (see 
Fig. [3]). This agreement unambiguously demonstrates 
that our treatment of non-universal corrections is reli- 
able in the simulated parameter range (linear fits for 
C < 0.75). Away from unitarity we find T c /e F = 0.202(9) 
and 0.252(15) for k = 0.217(2) and 0.474(8), respectively. 

The results for the strongly correlated regime essen- 
tially determine the general shape of the universal curve 
[T c /ef](k) shown in Fig. Q] Deep in the BEC regime 
(k ^> 1) the critical temperature is that of a weakly in- 
teracting Bose gas of strongly bound dimers which is ex- 
pected to increase on approach to the resonance. In the 
BCS limit (k < 0, 3> 1) the T c -curve starts from 
exponentially small values for k — > — oo, and thus the 
crossover between the two limiting regimes necessarily 
features a maximum in T c /ef- The results in Fig. [T] 
clearly show that this maximum must be on the BEC 
side (k > 0). The value at the maximum appears to be 
surprisingly high. For comparison, we show in Fig. [1] the 
critical temperatures of a Bose gas with hard- and soft- 
core sphere potentials with scattering length as = 0.6a. 
The tremendous computational cost required to deter- 
mine each point in the crossover regime reliably did not 
allow us to precisely locate the position of the maximum 
in the kpa ~ 1 region. 

The behavior of the critical temperature on the BEC 
side revealed by our simulations, suggests that the short- 
range structure of the strongly correlated state is radi- 
cally different from that of the compact-molecule Bose 
gas (and, obviously, also from that of the BCS state) in 
a broad range of n. In other words, we are dealing with 
two crossovers — one is from BCS to the substantial uni- 
tarity regime and the other is from the unitarity regime 
to BEC. 

To summarize, we performed first-principle simulations 
of the two-component unpolarized Fermi gas with res- 
onant inter-particle interaction obtaining the universal 
critical temperature T c jep = 0.152(7) at the unitarity 
point k = l/fci?a = thereby resolving the earlier con- 
troversy between the results of Refs. [15] and [lj|. We 
also obtain T c away from unitarity on the BEC side al- 
lowing one to sketch the general dependance [T c /ep](K) 
with a maximum on the BEC side, in a good quantitative 
agreement with the mean-ficld-type prediction by Hauss- 
mann et al. 0]. After our results were announced 21 1, 
the Seattle group reconsidered their previous estimate of 
the critical temperature 22]. The new results are in ex- 
cellent agreement with the values claimed here both at 
and away from unitarity. 

The simulations were performed on the supercomput- 
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